library(tidyverse)
library(here)
library(plotly)
library(htmlwidgets)
load(file = here("data/data_tidy.RData")) # load object data_tidy
BodyLengths <- tibble(Species = c("Ceriodaphnia sp.", "C. megalops", "D. cucullata", "D. curvirostris", "D. hyalina var. gellata", "D. hyalina var. lacustris", "D. hyalina", "D. longispina", "D. pulex", "D. magna", "Calanoid copepods", "Cyclops", "Bosmina coregoni", "B. longirostris", "Sida sp.", "S. crystallina", "Chydorus ovalis", "Eurycercus lamellatus", "Alona spp.", "A. quadrangularis", "Asplanchna", "Keratella spp.", "Ostracod", "Diaptomus"),
BodyLength = c(0.4636, 0.925, 0.8593, 1.58, 0.9064, 0.9064, 0.9064, 1.0347, 1.1656, 1.4214, 0.6860, 1.0369, 0.4230, 0.3637, 0.5075, 0.5075, 0.475, 1.975, 0.5235, 0.9679, 0.4747, 0.3495, 1.25, 0.9886)) # species and their body sizes
data <- left_join(data_tidy, BodyLengths, by = "Species") # add species body size to the dataframe
data <- mutate(data, SizeClass = case_when(BodyLength <= 0.6 ~ "Small (<= 0.6 mm)",
BodyLength <= 1.0 ~ "Medium ( 0.6 < x <= 1.0 mm)",
BodyLength > 1.0 ~ "Large (> 1.0 mm)"))
plotlist <- list()
lakes <- unique(data$LakeName)
for (i in 1:length(lakes)){
Ldata <- filter(data, LakeName == lakes[i])
if (nrow(Ldata) > 0) {
summary <- Ldata %>% group_by(Species, SizeClass, Date) %>% summarize(sum_counts = sum(Counts, na.rm = TRUE),
sum_total = sum(Total_individuals, na.rm = TRUE)) # Calculate counts per Date
summary <- summary %>% group_by(SizeClass, Date, sum_total) %>% summarize(sum_counts = sum(sum_counts, na.rm = TRUE))
summaryV <- Ldata %>% group_by(Species, Date) %>% summarize(sum_counts = sum(Counts, na.rm = TRUE),
sum_total = sum(Total_individuals, na.rm = TRUE))
summaryV <- summaryV %>% mutate("Relative_abundance" = sum_counts / sum_total * 100) # Calculate relative abundance per Date
summaryW <- left_join(summaryV, BodyLengths, by = "Species")
WeightedMean <- summaryW %>% group_by(Date) %>% summarize(WeightedMean_mm = weighted.mean(BodyLength, Relative_abundance))
#print(knitr::kable(WeightedMean,
#format = "markdown",
#caption = paste0("The weighted mean for ", lakes[i], " in ", years[j])))
summary <- summary %>% mutate("Relative_abundance" = sum_counts / sum_total * 100) # Calculate relative abundance per Date
scaleFactor <- max(summary$Relative_abundance) / max(WeightedMean$WeightedMean_mm)
result <- ggplot()+
geom_col(data = summary, aes(x = Date, y = Relative_abundance, fill = SizeClass), color = "black", linewidth = 0.05)+
geom_point(data = WeightedMean, aes(x = Date, y = WeightedMean_mm * scaleFactor), color = "darkmagenta")+
scale_y_continuous(name = "Relative abundance (%)", sec.axis = sec_axis(~./scaleFactor, name = "Weighted mean in mm"))+
theme_minimal()+
theme(
axis.title.y.right=element_text(color = "darkmagenta"),
axis.text.y.right=element_text(color = "darkmagenta"))+
labs(x = "Date", title = paste0("Relative abundance of each size class and the weighted means per date of ", lakes[[i]]))
plotlist[[(length(plotlist)+1)]] <- result # add the plot to the list of plots
} else {
warning(paste0("No data for ", lakes[[i]]))
}
}
htmltools::tagList(lapply(1:length(plotlist), function(x) { ggplotly(plotlist[[x]]) }))